Van der Waals Density Functional for Layered Structures 
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To understand sparse systems we must account for both strong local atom bonds and weak 
nonlocal van der Waals forces between atoms separated by empty space. A fully nonlocal functional 
form [H. Rydberg, B.I. Lundqvist, D.C. Langreth, and M. Dion, Phys. Rev. B 62, 6997 (2000)] 
of density-functional theory (DFT) is applied here to the layered systems graphite, boron nitride, 
and molybdenum sulfide to compute bond lengths, binding energies, and compressibilities. These 
key examples show that the DFT with the generalized-gradient approximation does not apply for 
calculating properties of sparse matter, while use of the fully nonlocal version appears to be one 
way to proceed. 
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Calculations of structure and other properties of sparse 
systems must account for both strong local atom bonds 
and weak nonlocal van der Waals (vdW) forces between 
atoms separated by empty space. Present methods are 
unable to describe the true interactions of sparse systems, 
abundant among materials and molecules. Key systems, 
like graphite, BN, and M0S2, have layered structures. 
While today's standard tool, density-functional theory 
(DFT), has broad application, the common local (LDA) 
and semilocal density approximations (GGA) []]■ Q. 0. |~| 
for exchange and correlation, E^ c [n], fail to describe the 
interactions at sparse electron densities. Here we show 
that the recently proposed density functional [j| with 
nonlocal correlations, Ejr[n], gives separations, binding 
energies, and compressibilities of these layered systems in 
fair agreement with experiment. This planar case bears 
on the development of vdW density functionals for gen- 
eral geometries 0, Q , as do asymptotic vdW functionals 

a. 

Figure 1 with its 'inner surfaces' defines the problem: 
voids of ultra-low density, across which electrodynamics 
leads to vdW coupling. This coupling depends on the po- 
larization properties of the layers themselves, and not on 
small regions of density overlap between the layers, ex- 
cluding proper account in LDA or GGA. For large inter- 
planar separation d the vdW interaction energy between 
planes behaves as —C4/d 4 , while LDA or GGA necessarily 
predicts an exponential falloff. Layers rolled up to form 
two (i) nanotubes with parallel axes a distance I apart 
interact as — C5/Z 5 , or (ii) balls (e.g., Ceo), a distance r 
apart, as —cg/r 6 . If by fluke an LDA or GGA were to 
give the correct equilibrium for one shape, it would neces- 
sarily fail for the others. The simple expedient of adding 
the standard asymptotic vdW energies as corrections to 
the correlation energy of LDA or GGA also fails. The 
true vdW interaction between two close sheets must be 
(i) substantially stronger (Fig. 1), (ii) seamless, and (iii) 



saturate as d shrinks (Fig. 2). 

Like earlier work directly calculating nonlocal corre- 
lations between two jellium slabs [jj], the vdW density 
functional (vdW-DF) theory || used here exploits as- 
sumed planar symmetry. It divides the correlation en- 
ergy functional into two pieces, E c [n] = E® [n] + [n] , 
where -E^M ^ s defined to include the longest ranged or 
most nonlocal terms that give the vdW interaction and 
to approach zero in the limit of a slowly varying den- 
sity. The term E®[n] is also nonlocal, but approaches 
the LDA in this limit. The two terms are approximated 
differently. Here the LDA is used for E® [n], as the LDA 
should be much more accurate with the principal longest 
range terms separated off. Ultimately we plan to derive 
a gradient or GGA expansion appropriate to £^?[n], but 
the resulting corrections are expected to be small. Thus 
we use 



E c [n]*E^ A [n] 



E?\n] 



(1) 



Long range terms are less sensitive to details of the sys- 
tem's dielectric response. Thus very simple approxima- 
tions for the dielectric function are made for [n]. We 
care to make the polarization properties of a single layer 
come out reasonably. 

For systems with planar symmetry, the predominant 
component of the nonlocal correlation energy Z?" 1 giving 
the vdW forces can be determined simply by comparing 
the solutions of the Poisson equation V • (eV$) = and 
the Laplace equation V 2 <!> = |5(. The crux of the ap- 
proximation is the use of a simple plasmon-pole model 
for the dielectric function, 



ek(z, iu) = 1 + 



u* + {v F (z)q k f/Z + qi/A' 



(2) 



where w 2 (z) = Airn(z)e 2 /m, and mvp(z) = (37r 2 n(z)) 1 / 3 
are functions of the local density n(z). The 2D wave 
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FIG. 1: Schematics of the vdW forces in sparse, layered mate- 
rials. Top panel: Laterally averaged electron-density profile of 
two graphene layers at equilibrium separation d ~ 3.5 A. The 
dotted ellipse indicates the region determining the interlayer 
interaction within LDA or GGA. It is unable to describe the 
true long range interaction (horizontal arrow) between corre- 
lated charge fluctuations (extrema illustrated by those of the 
finite-frequency charge response, and labelled by vertical ar- 
rows) . Bottom panel: Full (solid curve) vs. the asymptotic 
form of -B" 1 — * —C4/0I 4 (dashed line). Also shown: full ground 
state interaction energy (dash-dotted line). Insets: Sketches 
of the correlated charge fluctuations for different layer separa- 
tions, explaining why l-E" 1 ! > 04 /d 4 . In the asymptotic region 
(top inset), attractive interactions BC and AD are nearly can- 
celled by the repulsive interactions AC and BD, leaving little 
more than the weak dipolar interaction. Near equilibrium 
separation (bottom inset), the repulsive interactions AC and 
BD remain specified by the layer separation d, but the rela- 
tive strengthening of attractive interaction BC substantially 
overcompensates the relative weakening of the attractive in- 
teraction AD. At still smaller d, i?" 1 becomes weaker and 
saturates (Fig. 2, inset). 



vector perpendicular to the z direction is k, iu is the 
imaginary frequency and q\ = k 2 + q 2 ^_ mimicks the 3D 
wave vector, with the z direction accounted for by the 
n(z) dependence and q± taken as a constant. 

Our scheme, applied to jellium surfaces in Rcf. |5|, is 
used here for layered systems, for simplicity a two-layered 
system (e.g., two parallel graphene sheets) arranged per- 
pendicular to the z axis, at positions 2 = and z = d, 
respectively. The first step is to calculate the density in 
DFT with a suitable approximation, which we take to be 
an appropriate flavor of the GGA, and to average it in 
the lateral direction to yield n(z). The next two steps 
represent the key to the success of our approximation. 
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FIG. 2: Calculated vdW-DF results for the interlayer binding 
(adhesion) between two sheets of graphene (upper panel) and 
M0S2 (lower panel) as functions of the interlayer separation d. 
In each panel it is compared (solid curves) with the underlying 
GGA (dashed curves) , and with the adhesion in the absence of 
the vdW-DF term E™ 1 (dash-dotted curves) calculated using 
only the first two terms of Eq. @ for E xc . The nonlocal- 
correlation vdW contribution E™ 1 alone is shown separately 
(dotted curves). The figure shows that no binding results from 
exchange or local correlation; the vdW interactions provides 
adhesion in fair agreement with experiment (Table I). The 
inset shows that our calculated value of E^ 1 saturates for small 
d values. 

First one calculates from first principles the perpendicu- 
lar static polarizability a of a single layer of the material 
(a graphene layer in this example) in the ground-state 
DFT scheme with the above GGA flavor. Then one fixes 
the constant q± so that dielectric model @ gives pre- 
cisely the same polarizability, that is 

1 f \ 1 1 

a = — dz 1 — — . (3) 

Air J [ e (2,0)_ 

This step mitigates the errors introduced both by the 
lateral averaging and the dielectric model. The leading 
vdW interaction is proportional to an integral over the 
square of the dynamic polarizabilities. Previous work 
in the asymptotic limit |a| shows that simple, properly 
scaled dielectric functions giving the correct static po- 
larizabilities form good approximations to the dynamic 
ones, hence giving accurate vdW interactions. 

The correlation energy can be calculated from the 
charge response |l0j . With planar symmetry only the 
response to an arbitrarily charged sheet p oc exp(zk ■ 
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FIG. 3: Comparison of the vdW-DF layer-to-layer binding 
in the staggered bulk hexagonal-BN system (solid lines) and 
adhesion between two sheets of BN (broken lines) shown as 
functions of the interlayer separation d — c/2. In both de- 
scriptions the GGA gives unphysical results for bond and 
equilibrium separation. The total adhesion is dominated by 
the contribution of the vdW interaction . The nearest- 
neighbor interlayer binding (broken lines) is seen to strongly 
dominate the bulk adhesion (solid lines) for layered materials. 



r)S(z — Z-) placed at an offset z_ (z_ <C 0) from one end 
of the sample needs to be calculated. Integrating Pois- 
son's equation with e as in @ one finds the z component 



of the electric field E z (z + ) at z 
given by || 



» d. Then £?' is 



E 



ni 



-A 



'du 
2^ 



d 2 k E z (z+) 
In 



(2nf E°(z + Y 



(4) 



where A is the lateral area and E®(z + ) is the electric 
field component from the charge sheet at z — z_ in the 
absence of the sample. In Eq. Q the spatial dependences 
of the two E's cancel in the ratio, leaving only the k 
dependence. For each of the bulk solids considered here 
we apply this method to 32 and 30 layer slabs to give, by 
subtraction, a well converged per layer value of £7" . 

The approximation for the exchange energy functional 
-Exjn] should be consistent with our approximation for 
E®[n], that is a local or semilocal one. We use a ver- 
sion that best approximates exact density functional ex- 
change, namely Zhang and Yang's 0] (ZY) "revPBE" ex- 
change functional E^ Y [n], where the parameter control- 
ling the large gradient limit in the PBE exchange func- 
tional is fitted to exact density-functional exchange 
calculations. Unlike exchange m PW91 0], PBE 0, 
and RPBE E^ Y shows no tendency to bind any of 
the vdW systems we have tried it on, unless the vdW 
correlation is actually included in the calculation. We 
thus take 



E xc [n] = E? [n] + E^ UA [n] + E» l [n] 



(5) 



The first two terms form an important and recommended 
starting point for testing nonlocal vdW functionals. They 



TABLE I: Properties of graphite, BN, and M0S2, calculated 
with vdW-DF and compared to experimental data in paren- 
theses, when available. The table shows the geometry (a, c), 
binding energy (-En), bulk modulus (-Bo), and elastic constant 
(C33). GGA (not shown) binds very weakly at unphysically 
large c or not at all, depending on substance and GGA flavor. 



Graphite 



BN 



M0S2 



a [A] 
c[A] 

E [meV/at] 
B [GPa] 
C33 [GPa] 



2.47 (2.46)" 
7.52 (6.70) d 
24 (35 ± 10) e 

12 (~33) a 

13 (37-41)" 



2.51 (2.50) 6 
7.26 (6.66)' 

26 

11 

11 



3.23 (3.160) c 
12.6 (12.29) c 

60 

39 

49 



a Ref. 
6 Ref. 
c Ref. 



d Ref. 2£ 
E Ref. 27 



assure a vdW attraction that actually comes from terms 
that in principle are capable of giving it. 

The GGA is used in order to have a good account of 
strong valence bonds and densities of individual layers, 
and the above scheme is implemented by sclf-consistently 
calculating the energy E GGA and density n GGA , typically 
in the revPBE flavor of GGA. The total energy, as a func- 
tion of the lattice parameters, is then calculated as E 

£GGAr n GGA] + /^^GGA^ where A [ n ] = £nl[ n ] _ E[ 



eradr 



and£s rad f 



E GGA [r 



E, 



LDAr 



This approximation 



treats A as a post GGA perturbation. The GGA calcula- 
tions are done using the plane-wave pseudopotential [l2T | 
DFT code dacapo O. 

The three materials considered here have layered struc- 
tures with a soft direction perpendicular to the lay- 
ers 01- Strong covalent bonds occur between the 
atoms within the hexagonal layers (lattice constant a), 
while weak vdW interactions occur between the layers. 
Graphite has a staggered or A-B type stacking with the 
intralayer C atoms in regular, planar hexagons stacked 
with corners on top of centers. Boron nitride, isoelec- 
tronic with graphite, has no shift in hexagon locations 
but B and N atoms in alternate positions along the re- 
direction. A similar alternation occurs between Mo and 
S atoms in the A-B type structure of M0S2 0] • 

Calculations are also performed for bilayers of the 
three systems. The binding-energy curve of two paral- 
lel graphene sheets 0], i.e. the total-energy difference 
at separation d and at infinite separation, respectively, 
for varying d (Fig. 2) gives in the approximation [f| with 
vdW-DF, Eq. JSJ), a close relation to experimental find- 
ings for binding energy and equilibrium distance. The 
GGA curve, on the other hand, is completely wrong, 
which we blame on absence of vdW effects. The inset 
shows that our E^r contribution gives stronger binding 
than the traditional asymptotic vdW interaction. The 
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binding-energy curve for two parallel BN sheets (Fig. 3) 
differs little from that of bulk BN as a function of lattice 
constant c. At equilibrium the calculated E^ 1 difference 
is ~ 4%, attributable mostly to 2nd nearest layers in the 
bulk material 17]. 

Bulk-modulus values are computed to geth er with other 
structure (a,c) and bonding properties [18j (Table I). It 
is crucial to densely sample the region of (a,c) values 
around the optimal structure (do, Co), and we use a new 
method 0] for direct evaluation of both structure, bind- 
ing energy, and bulk modulus Bq in the relevant range 
of a and c values [i(J. GGA values were also computed 
but not shown. For all three materials the various GGA 
flavors give no binding or bind very weakly at unreason- 
ably large separation. The vdW-DF, on the other hand, 
gives values for lattice parameters and cohesive energy in 
good agreement and for bulk modulus in fair agreement 
with experimental values, when available. 

In conclusion, we recommend the replacement of GGA 
as a standard method in total-energy calculations with 
vdW-DF as given in Eq. for descriptions of layered 
systems that contain sparse electron distributions. This 
will give the right qualitative character of soft bonds, 
including saturation of the vdW potential at small sepa- 
rations, and even quantitative ones, like physical values 
of bond lengths, binding energies, and compressibilities 
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